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ABSTRACT 

Context. It is a well established fact that some YSO jets (e.g. RW Aur) display different propagation speeds between their blue and 
red shifted parts, a feature possibly associated with the central engine or the environment in which the jet propagates. 
Aims. In order to understand the origin of asymmetric YSO jet velocities, we investigate the efficiency of two candidate mechanisms, 
one based on the intrinsic properties of the system and one based on the role of the external medium. In particular, a parallel or anti- 
parallel configuration between the protostellar magnetosphere and the disk magnetic field is considered and the resulting dynamics 
are examined both in an ideal and a resistive magneto-hydrodynamical (MHD) regime. Moreover, we explore the effects of a potential 
difference in the pressure of the environment, as a consequence of the non-uniform density distribution of molecular clouds. 
Methods. Ideal and resistive axisymmetric numerical simulations are carried out for a variety of models, all of which are based 
on a combination of two analytical solutions, a disk wind and a stellar outflow. The initial two-component jet is modified either by 
inverting the orientation of its inner magnetic field or imposing a constant surrounding pressure. The velocity profiles are studied 
assuming steady flows as well as when strong time variable ejection is incorporated. 

Results. Discrepancies between the speeds of the two oppositely directed outflows can indeed occur both due to unaligned magnetic 
fields and different outer pressures. In the former case, the asymmetry appears only on the dependence of the velocity on the cylindrical 
distance, but the implied observed value is significantly altered when the density distribution is also taken into account. On the other 
hand, a non-uniform medium collimates the two jets unevenly, directly affecting their propagation speed. A further interesting feature 
of the pressure-confined outflow simulations is the formation of static knots whose spacing seems to be associated with the ambient 
pressure. 

Conclusions. Jet velocity asymmetries are anticipated both when multipolar magnetic moments are present in the star-disk system 
as well as when non-uniform environments are considered. The latter case is an external mechanism that can easily explain the large 
time scale of the phenomenon, whereas the former one naturally relates it to the YSO intrinsic properties. 

Key words. ISM: jets and outflows - MHD - Stars: winds, outflows, pre-main sequence 



1. Introduction 

Over the last few years, the two-component jet scenario emerges 
as a strong candidate for describing Young Stellar Object (YSO) 
outflows. Obser vation al data of Cl assical T Tauri Stars (CTTS) 
(Edwards et al. 2006; Kwan et al. 2007 ) indicate the presence 
of two genres of winds: one being ejected radially out of the 
central object (e.g. Sauty & Tsinganos |1994t Trussoni et al. 
I1997t Matt & Pudritz 12005b and the other being launched at 
a constant angle with respect to the disk plane (e.g. Blandford 
& Payne [T9821 Tzeferacos et al. [20091 Salmeron et al. [20TTT) . 
Consequently, CTTS outflows may be associated with either a 
stellar or disk origin, or with both outflow components having 
comparable contributions. In addition, such a scenario is sup- 
ported by theoretical arguments (Ferreira et al. 12006b . An ex- 
tended disk wind is required for the explanation of the observed 
YSO mass loss rates, whereas a pressure driven stellar outflow 
is expected to propagate in the central region, possibly account- 
ing for the protostellar spin down (Matt & Pudritz 2008; Sauty 



et al. 1201 lb . Numerical simulations have been also employed re- 
cently to investigate the various aspects of two-component jets 
(Meliani et al. 125551 Fendt[2559l Matsakos et al. [25591 hereafter 
M09). 

Detailed observations of some YSO with bipolar flows have 
shown peculiar velocity asymmetries between the blue and red 
shifted regions (e.g. Woitas et al. 2002; Coffey et al. 2004; Perrin 
et al. 12007b . In particular, recent estimates of the RW Aur jet 
indicate that the speed of the approaching lobe is roughly 50% 
higher as compared to the receding one (Hartigan & Hillenbrand 
2009). Although the above studies suggest that the asymmetry 
originates close to the source, observations of the same object 
carried out by Melnikov et al. (2009) point out a similar mass 
outflow rate for the two opposite jets, concluding that their dif- 
ferent speeds are due to environmental effects. Thus, it is still an 
open question whether the discrepancy between the two hemi- 
spheres relies on an intrinsic property, such as the magnetic field 
configuration, or an external factor, such as the physical condi- 
tions of the surrounding medium. 
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Spectropolarimetric measurements of T Tauri stars suggest 
multipolar magnetospheres, with the dipolar component n ot al- 
ways parallel to the rotation axis (Valenti & Johns-Krull 2004; 
Daou et al. 20061 Yang et al. 125571 Donati & Landstreet l2559l) . 
Non-equatorially symmetric field topologies may affect both the 
way that matter accretes on the protostar (Long et al. 120071 
Long et al. 12008b and the jet launching mechanisms. For jets, 
it is not clear whether the magnetic field asymmetry persists 
for timescales comparable to the jet propagation timescales 
(decades-centuries). In fact, a periodic stellar activity could 
be directly associated with jet variability. In a similar context, 
Lovelace et al. (2010) simulated the central region of YSOs as- 
suming complex magnetic fields. In one interesting case, where 
the magnetosphere is a combination of dipolar and quadrupo- 
lar moments, they find the extreme scenario of one-sided con- 
ical outflows being ejected from the star-disk interaction inter- 
face. On the other hand, the accretion disk could also have a 
quadrupolar magnetic field originating either from the dynamo 
mechanism or from advection (Aburihan et al. 12001b . In fact, it 
has been shown that non-bipolar disk field topologies can be a 
potential source for asymmetry in AGN jets (Wang et al. 119921 
Chagelishvili et al. [T996b . 

The external medium may play a relevant role in the dynam- 
ical propagation of the outflow; we refer again to the detailed 
HST observations of RW Aur by Melnikov et al. (120091) , whose 
findings can be consistent with the effects of inhomogeneities 
in the environmental conditions. Clumpy and filamentary struc- 
tures are observed in molecular clouds over a wide range of 
scalelenghts, implying the presence of an external medium that 
surrounds the jets. In addition, an environment affecting jet 
propagation may be present but remain almost undetectable in 
observations, as discussed recently by Te§ileanu et al. (120121) . 
Theoretical arguments suggest that density anisotropies in the 
vicinity of YSOs could collimate YSO outflows and even induce 
oscillations in the jet's cross section (Konigl 1982). Numerical 
simulations have investigated the effects of the surrounding gas, 
such as a collapsing environment whose ram pressure colli- 
mates a spherical wind (Delamarter et al. 2000) or an isother- 
mal medium with a toroidal density distribution that results in 
a cylindrically shaped outflow (Frank & Noriega-Crespo [T9941 
Frank & Mellema 1996 ). Another class of simulations has stud- 
ied the collimating role of a vertical outer magnetic field, whose 
magnetic pressure effectively confines an isotropic stellar wind 
(Matt & Bohm l2053l Matt et al. [2003b . We note that although 
the external pressure is not thought to be responsible for the jet 
collimation observed at larger scales, it could still affect it and 
hence modify the propagation speed. 

The goal of the present work is to study the feature of the 
velocity asymmetry within the context of the two-component jet 
models presented in M09. We take advantage of both theoretical 
and numerical approaches (Gracia et al. 120061) , setting as initial 
conditions a combination of two analytical YSO outflow solu- 
tions. Tabulated data of such MHD jets have been derived using 
the self- similarity assumption (Vlahakis & Tsinganos 1998 ) and 
are available for both disk winds and stellar jets. We employ 
the Analytical Disk Outflow (ADO) solution from Vlahakis et 
al. (120001) a nd the Analytical Stellar Outflow (ASO) model from 
Sautv et al. (120521) . 

Matsakos et al. (2008| have addressed the topological sta- 
bility, as well as several physical and numerical properties of 
each class of solutions separately. In M09 the two comple- 
mentary outflows were properly mixed inside the computational 
box, such that the stellar outflow dominated the inner regions 
and the disk wind the outer. The stability and co-evolution of 
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Fig. 2. The northern hemisphere assumes a lower pressure, and 
hence a lower degree of collimation, as compared to the south- 
ern. The sketch is not in scale, the effects of this mechanism 
will hold at any distance wherein distinct environments are en- 
countered. The magnetic field topology here is simple, a dipolar 
stellar magnetosphere is surrounded by a bipolar disk field. 



several dual component jet cases was investigated as a func- 
tion of the mixing parameters and the enforced time variabil- 
ity. Furthermore, the introduction of flow fluctuations generated 
shocks, whose large scale structure demonstrated a strong re- 
semblance to real YSO jet knots. 

In this paper, two main scenarios are examined as possi- 
ble mechanisms to produce velocity asymmetries. The first one 
refers to intrinsic YSO properties introducing distinct magnetic 
field topologies at each side of the equator. The second one refers 
to external effects arising from differences in the ambient pres- 
sure. 

In the first case, presented in Fig.[T] each hemisphere is con- 
sidered to have one of the following two magnetic field configu- 
rations: the field lines coming out of the stellar surface are par- 
allel (north) or anti-parallel (south) to the large scale magnetic 
field. This implies a quadrupolar disk field and a dipolar magne- 
tosphere (left), or equivalently a bipolar disk field surrounding a 
magnetic quadrupole (right). 

In the second scenario, shown in Fig. [2) the properties of 
the medium along the propagation axis are assumed to affect the 
collimation of the jet. Specifically, the outer pressure could mod- 
ify the diameter of the outflow and in turn it would determine 
the wind speed. This mechanism does not depend on any YSO 
time scale, and it could appear at any height that the jet encoun- 
ters heterogeneous environment. We note that although a thermal 
pressure of a hot tenuous medium is used in the simulations, this 
could equivalently represent the pressure due to a cold dense en- 
vironment, turbulence or even a large scale magnetic field. Stute 
et al. (120081) truncated ADO solutions imposing a weaker out- 
flow of the same type at outer radii. Here, we extend their study 
assuming a constant surrounding pressure instead, also with dif- 
ferent values at each jet direction. 

All our models focus on large enough scales avoiding the 
complicated dynamics of the star-disk interaction (e.g. Bessolaz 
et al. 120081 Zanni et al. l2059l Lovelace et al. l20T0b . Essentially, 
these mechanisms are included merely as boundary conditions, 
an approach frequently adopted in the literature of jet propaga- 
tion studies (e.g. Ouyed et al. 120031 Fendt 2006; Te§ileanu et al. 
125591) . 
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Fig. 1. Left panel: a quadrupolar disk field and a dipolar stellar magnetosphere. Right panel: a bipolar disk field and a quadrupolar 
stellar magnetosphere. The northern hemisphere of either case is referred to as a parallel configuration, whereas the southern as 
an anti-parallel. The dashed line indicates the surface that separates the disk wind from the stellar outflow. The simulations do not 
include the central region of the YSO, and hence we do not attempt to describe the star-disk interaction that takes place inside the 
grey ellipse. In addition, the requirement of V • B - means that field lines must return along the accretion disk. 



Simulations are carried out to study the time evolution as 
well as the potential final steady states, comparing the results 
obtained between the two sides of the equatorial plane. We are 
especially looking for emerging asymmetries in the vertical ve- 
locity profiles. The resistive MHD regime is of particular inter- 
est (Cemefjic et al. 2008) since reconnection plays a key role, 
mainly at the locations where the magnetic field inverts sign. 
Furthermore, all cases are re-examined assuming that the flow 
has a strong time variable character, an effect that produces dis- 
continuities along its axis. Note that we do not attempt to model 
the observed properties of the RW Aur jet, but rather to investi- 
gate in a general manner the intrinsic or extrinsic physical pro- 
cesses that could underlie the velocity asymmetry phenomenon. 

The paper is structured as follows. Section §[2 describes the 
jet models and provides information on the numerical setup. 
Section ^presents and discusses the obtained results while sec- 
tion §|4l reports the conclusions of this work. 



(47r)" 1/2 , satisfies the condition V • B = 0. Moreover, the electric 
current is defined as / = V x B, and has also absorbed a factor 
(47r)" 1/2 . The gravitational potential is given by O = -GM/R, 
where G is the gravitational constant, M the mass of the cen- 
tral object and R the spherical radius. The magnetic resistivity 
tensor has been assumed diagonal, with all its elements equal to 
the scalar quantity 77. Finally, A represents the volumetric energy 
gain/loss terms and T = 5/3 is the ratio of the specific heats. 

We adopt the axisymmetric cylindrical coordinates (r, z) and 
we will express all quantities in code units. A code variable U 
is converted to its physical units U' with the help of Uq, i.e. 
U' - UqU. The constants Uq are specified such that the values of 
the quantities in the final configuration of the simulations match 
the typical values of real jets, such as observed velocities and 
densities. Specifically, r$ = zo = 1.0 AU,po = 2.1 • 10" 18 gem -3 , 
V = 4.0 km s" 1 , P = 3.3 • 10" 7 dyne cm" 2 , B = 2.0 • 10" 3 G 
and to = 1.2 yr. 



2. Numerical models and setup 

2.1. MHD equations 

The non-dimensional resistive MHD equations written for the 
primitive variables are: 



g + V.(pV) = 0, 

^ + (V • V)V + -B x J + -VP = -VO , 

dt p p 

^ + V • VP + TPV • V = (r - 1)?]J 2 + A , 
dt 

dB 

— + V x (B x V) = -//V x / . 

dt 



(1) 

(2) 

(3) 
(4) 



The quantities p, P and V are the density, pressure and veloc- 
ity respectively. The magnetic field B, which includes the factor 



2.2. Analytical solutions 

The employed ADO solution, denoted with subscript D, is ap- 
propriate for describing disk winds. It is a steady state radially 
self- similar MHD outflow, namely, each quantity has a certain 
scaling along the disk radius based on the Keplerian rotation 
profile. Therefore, if all flow variables are known along a par- 
ticular fieldline the spatial distribution of the entire solution can 
be reconstructed. This type of models have conical critical sur- 
faces and consequently diverge on the axis. However, we sub- 
stitute the inner region with the ASO solution, denoted with S, 
a meridionally self-similar jet. In this case, the critical surfaces 
have a spherical shape that effectively models stellar outflows. 
Some kind of energy input is required at the base of the wind 
to accelerate the plasma, but we do not treat this region because 
we simulate high altitudes that correspond to its super Alfvenic 
domain. Since the ADO solution is described by the poly tropic 
index T - 1.05, we specify the same value in the simulations 
which implies an energy source term A = (5/3 - T)PV • V. We 
have verified that this does not significantly affect the dynamics 
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(see also Matsakos et al. 2008 ). In other words, although ther- 
mal effects are thought to be crucial for the mass loading close 
to the disk surface (Ferreira & Casse 2004 ), we can neglect them 
above it where the magneto-centrifugal mechanism takes over. 
The adopted analytical solutions together with their parameters 
are presented in detail in M09. 

2.3. Normalization and mixing 

The two-component jet model parameters can be classified into 
two categories. The first one describes the relative normalization 
of the analytical solutions with the help of the following ratios: 



R* 



Vs* 
Vd* 



(5) 



where R* is the Alfvenic spherical radius of the ASO model, r* 
the cylindrical radius of the Alfvenic surface on a specific field- 
line of the ADO solution and the star denotes the characteris- 
tic values at these locations. The subscripts L, V and B stand 
for length, velocity and magnetic field, respectively. We specify 
iy - 5.96, such that the two solutions correspond to the same 
protostellar gravitational field. Following the guidelines of M09 
we set = 0.1. Note that due to the non-existence of a charac- 
teristic length scale in the ADO solution, the value of this param- 
eter does not play any role in the combination of the two wind 
components. We also define £^ =0.01, in order to provide com- 
parable magnetic fields for the two outflows in the initial con- 
ditions of the numerical setup. Larger or lower values for this 
parameter would result in either a dominant stellar outflow and 
a negligible disk wind or the opposite, respectively. 

The second class of parameters controls the combination of 
the two analytical models. In order to achieve a smoothly varying 
magnetic field and ensure the dominance of the protostellar mag- 
neto sphere in the central region, we follow an improved mixing 
process as compared to the one presented in M09. The merging 
depends on the magnetic flux functions A D and As, which label 
the field lines of each initial one-component solution. Once a to- 
tal magnetic flux A is defined, then the poloidal component of 
the magnetic field is given by: 



£ p = (VAx0)/r, 



(6) 



where is the unit toroidal vector. However, a steep transition 
from one profile to the other could create artificial gradients in 
the flux which would be reflected in the magnetic field. As a 
result, B could end up being excessively strong or even inverted 
across the matching surface. 

Therefore, the footpoint of the fieldline that marks the tran- 
sition, A m i x , is chosen at a location where both Ad and As have 
roughly the same slope. Then, a first approximation of the to- 
tal magnetic flux is defined: A\ = A D + As + Ac, where Ac is 
a constant that helps normalize Ad and As. We require that the 
disk and stellar fields are exponentially damped in the central 
and outer regions respectively, which leads to an improved ap- 
proximation for the magnetic flux: 



1 - exp 



► A D + exp 



As. 



(7) 



Finally, with the help of A2 all two-component physical vari- 
ables, U, are initialized using the following gaussian-type mix- 
ing function: 



U 



, 1 " ex P 



A 2 



£/d + exp 



(8) 



This relation is also applied to generate the total initial two- 
component magnetic flux, A. The poloidal field is then derived 
from Eq. © and it is divergence free by definition. 

Moreover, since accretion and protostellar variability are ex- 
pected to introduce fluctuations in the ejection, the velocity pre- 
scribed on the bottom boundary is multiplied with the following 
function: 



/s(r,f) = 1 +/7sin(^)exp 

U var / 



r 

F var 



(9) 



where r var = 2.5 is the period of the pulsation and r var = 5 is 
roughly the cylindrical radius at which the matching surface in- 
tersects the lower boundary of the computational box. The frac- 
tional variability is set to p - 0.5 in the time variable models, 
in order to check the stability in the presence of strong pertur- 
bations and enhance reconnection when resistivity is included. 
Essentially, Eq. © provides a gaussian spatial distribution of a 
periodic velocity variability of ~3 yr, that produces shocks and 
knot-like structures. 

2.4. Numerical models 

The anti-parallel configuration is initialized assuming a sud- 
den inversion of the magnetic field beyond the fieldline rooted 
at ri nv = 5. In non-resistive MHD the two systems are en- 
tirely equivalent due to the invariance of the axisymmetric MHD 
equations under the following transformations: i) B — > —B, ii) 
— > -Bq and — > -V^. The finite values of the electric 
conductivity and its dependence on quantities such as J is an 
unknown factor in YSO outflows. Therefore, we follow a simple 
approach applying a constant physical resistivity, with the value 
77 = 1.0, an approximation adequate to provide a rough estimate 
of how the system evolves in the presence of non-ideal effects. 
As explained in this value can be considered small because it 
ensures the dominance of the advection terms over the resistive 
ones by several orders of magnitude throughout most parts of the 
computational box. However, this does not hold true on recon- 
nection separatrices and hence its impact cannot be neglected 
there. Finally, we note that even in the ideal MHD cases, the 
finite size of the grid cells inevitably introduces artificial recon- 
nection. We minimize its effects by increasing the resolution and 
we estimate its influence by comparing results of different mesh 
refinements. 

On the other hand, the truncation of the jet for the low pres- 
sure hemisphere is initialized with the following process. We cal- 
culate the total pressure of the outflow, thermal plus magnetic, at 
the location rp = 30 on the bottom boundar}Q. Then, a constant 
pressure P^ of the same magnitude is applied on the domain 
beyond the limiting fieldline rooted at rp. In addition, the den- 
sity is arbitrarily reduced by two orders of magnitude, whereas 
the velocity and magnetic fields are set to zero. Although pres- 
sure equilibrium holds at the bottom boundary, the weaker field 
and jet pressure at higher altitudes cannot compensate for the 
push of the external medium. As a result, a radial collapse is ex- 
pected during the first steps of the simulation until the horizon- 
tal forces are balanced. For the jet propagating in the opposite 
direction, we assume another simulation case where the trunca- 
tion is still at rp, but a higher outer pressure is applied with the 
value Pn = 2P^. Finally, we include a case where the medium 
has Pyn = 6/Y in order to discuss the morphology of the final 
steady state configuration. 



1 The plasma-/? (thermal over magnetic pressure) is approximately 
equal to 0.03 there. 
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Table 1. The numerical models, the name of which follows from 
the initials of their description. 



Name 


Variability 


MHD 


B orientation 


Outer P 


NIP 


Non-variable 


Ideal 


Parallel 


- 


NIA 


Non-variable 


Ideal 


Anti-parallel 


- 


NRP 


Non-variable 


Resistive 


Parallel 


- 


NRA 


Non-variable 


Resistive 


Anti-parallel 


- 


VIP 


Variable 


Ideal 


Parallel 


- 


VIA 


Variable 


Ideal 


Anti-parallel 




VRP 


Variable 


Resistive 


Parallel 




VRA 


Variable 


Resistive 


Anti-parallel 




NIPL 


Non-variable 


Ideal 


Parallel 


Low 


NIPH 


Non-variable 


Ideal 


Parallel 


High 


NIPVH 


Non-variable 


Ideal 


Parallel 


Very High 


VIPL 


Variable 


Ideal 


Parallel 


Low 


VIPH 


Variable 


Ideal 


Parallel 


High 


VIPVH 


Variable 


Ideal 


Parallel 


Very High 



Table \T\ lists the fourteen two-component jet models consid- 
ered, along with a brief description of their setup. In particular, 
models NIP and NIA focus on the effects of the parallel and 
anti-parallel magnetic field configuration in the ideal MHD for- 
mulation. Cases NRP and NRA include resistivity in order to 
understand the role of magnetic reconnection and Ohmic heat- 
ing, especially in the regions where the magnetic field inverts 
sign. On the other hand, different values of an external pressure 
are imposed in models NIPL, NIPH and VIPVH to investigate 
the effects of distinct environments. We also consider the time- 
variable version of all previous cases to examine their behavior 
in the presence of flow fluctuations, i.e. models VIP, VIA, VRP, 
VRA, VIPL, VIPH and VIPVH. 

2.5. Numerical setup 

The simulations are performed with PLUTC0, a modular 
shock-capturing numerical code for computational astrophysics 
(Mignone et al. 2007; Mignone et al.|2012). Cylindrical coordi- 
nates are chosen in 2.5D and second order accuracy is specified 
in both time and space. All simulations adopt the TVDLF solver 
(Lax-Friedrichs scheme), a choice that is essential to retain sta- 
bility but comes at the cost of increased numerical diffusivity. 
Outflow conditions are specified at the top and right boundaries, 
axisymmetric at the left, and all variables are kept fixed to their 
initial values on the bottom boundary. The size of our domain 
is < r < 100 in the radial direction and 10 < z < 310 in 
the vertical. We are interested in radial distances less that 50, 
but a larger box ensures that any boundary effects coming from 
the imposed zero derivative of B will not influence our results 
(Zanni et al. 120071) . Note that the detailed launching mechanisms 
of each component are not taken into account since the computa- 
tional box starts from a finite height. The grid resolution is set to 
500x 1500 and the simulations are carried out until t stop = 30. We 
point out that the non- variable and non-pressure-confined simu- 
lations reach an exact steady state before t = 10. In M09, we 
have studied the steady state of such two-component jet models 
in larger time scales and we have found that the configuration is 
time-independent to very high accuracy. The pressure-confined 
models require an order of magnitude larger evolution time to 
attain a quasi steady outcome, and hence we specify 4top - 300 
and a resolution of 300 x 900 for efficiency purposes. 



2 Freely available at http : //plutocode . ph . unito . it 
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Fig. 3. Logarithmic density contours along with magnetic field 
lines (left) and vertical velocity contours together with poloidal 
current lines (right) for the initial conditions of model NIP. The 
dashed line indicates the surface A m [ x around which the mix- 
ing is assumed, and also marks the polarity reversal for the anti- 
parallel configurations. The thick solid line shows the separa- 
trix beyond which a constant pressure is initially applied on the 
pressure-confined models. The initial setup of NIA is identical 
except for the opposite sign of B and J in the central regions. 



3. Results 

As an overall remark, all non-truncated models reach steady 
states in an equivalent physical time of a few years. The final 
configuration of the disk wind remains close to the initial setup, 
demonstrating the stability of the ADO solution despite its sig- 
nificant modification. On the other hand, the hot stellar com- 
ponent gets compressed around the axis, which also results in 
an increase of the flow speed there. Even in the cases where 
the flow is strongly fluctuating, the two-component jet retains 
on average the properties of the unperturbed model. A steady 
shock that forms along a somewhat diagonal line at the lower 
part of the box is found to causally disconnect the jet launching 
regions from the propagation domain. Although a different mix- 
ing has been assumed, all the above features are common to the 
two-component jet models presented in M09, wherein a detailed 
discussion can be found. 

3.1. Parallel/anti-parallel configurations in ideal MHD 

The left panel of Fig. [3] presents the density and the mag- 
netic field lines of the starting setup of NIP, whereas the right 
panel shows the distribution of the vertical velocity and the 
poloidal current contours. The initial configuration of NIA has 
an identical morphology apart from the reversal of B and / that 
takes place inside the region marked with a dashed fieldline. 
Essentially, this surface separates the hot stellar outflow from 
the cold magneto-centrifugal disk wind, that is supposed to carry 
most of the mass and angular momentum extracted from the sys- 
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Fig. 4. In the left panel, the image shows the logarithmic density 
and the magnetic field lines for the final steady states reached 
for models NIP (left hand side; parallel configuration) and NIA 
(right hand side; anti-parallel configuration). In the right panel, 
the vertical velocity distribution and electric current contours are 
shown for the same models. Each panel displays side by side 
both hemispheres of the YSO system, described in Fig. [TJ in 
order to compare the steady states reached by the two models. 
Matter accumulates along the current sheet of NIA reducing the 
density of inner radii as well. In addition, its velocity magnitude 
is similar to NIP although with a more extended radial profile. 



tern. Moreover, the pressure driven flow is faster as compared to 
the outer parts of the two-component model, a property associ- 
ated with the normalization and the acceleration efficiency of the 
specific analytical solutions employed. In relation to the discus- 
sion on the mixing function, Fig. [5] demonstrates that there is a 
restricted choice for the location of the matching surface due to 
the requirement that both solutions ought to have there a mag- 
netic field of a similar shape and magnitude. In other words, the 
roughly vertical field of the stellar component cannot be matched 
smoothly with the bent lines coming out of the disk at outer radii. 
Finally, the thick solid fieldline represents the truncation surface 
that is imposed in the simulations which study the effects of an 
external pressure. 

Each one of the two panels of Fig. [4] shows side by side the 
jet evolution of the parallel and anti-parallel configuration de- 
scribed in Fig. \T\ In particular, the left pair shows the density 
and magnetic field lines of the final steady states reached by NIP 
(left) and NIA (right). The flipping of the field is evident in the 
anti-parallel case, although the general structure is maintained. 
From the theoretical point of view of ideal MHD, the magnetic 
field reversal (B — > -B) of a flux tube would not change the dy- 
namics. However, the thickness of the current sheet that would 
form on its surface cannot be assumed zero in a physical situa- 
tion. 

The simulation of NIA gives a decreased value of density at 
inner radii and a peak along the field inversion surface. The latter 
feature forms simultaneously at all heights during the first steps 
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Fig. 5. The logarithmic distribution of the magnetic Reynolds 
number R m = rV z /r] (left panel) and Rp = (P/B 2 )R m (right panel) 
for both YSO hemispheres, i.e. NRP (on the left of each pair) and 
NRA (on the right of each pair). Note that their definition does 
not include the length scale of the magnetic field reversal, but we 
can still observe that the resistive effects are small over most of 
the computational domain. 



of the simulation, and it appears more prominent in lower reso- 
lutions. In fact, numerical reconnection inevitably appears when 
a curved current sheet is dragged through a square cell grid. The 
field, and in particular the strong toroidal component that domi- 
nates the magnetic pressure, are destroyed at the interface giving 
rise to a perpendicular force. Matter accumulates along that sur- 
face and eventually the thermal pressure compensates the lack 
of the magnetic one. More accurate solvers and higher resolu- 
tion would treat in a more consistent way the separatrix of the 
inversion of B. However, on the one hand the applied physical 
resistivity of the simulations presented below is well above the 
numerical diffusivity, and on the other, "smoothing out" effects 
are an important stability factor for the initial two-component 
magnetic field. 

The right pair of Fig. |4] plots the final vertical velocity dis- 
tribution and the poloidal currents for the corresponding models 
shown on the left. A wider radial profile of V z is observed in 
NIA as compared to NIP. Additionally, apart from the current 
sheet that forms in both models along the weak diagonal shock, 
another one also appears on the surface of the magnetic field 
reversal, as expected according to the previous discussion. 

3.2. Parallel/anti-parallel configurations in resistive MHD 

Based on the analysis of Cemeljic et al. (2008) we consider the 
magnetic Reynolds number R m = rV z /rj and the quantity Rp = 
(P/B 2 )R m , which measure the contribution of the resistive terms 
in the induction and energy equations, respectively. 

Fig. [5] displays log(7? m ) (left pair) and log(Rp) (right pair) 
for the final states of NRP (left) and NRA (right). Clearly, the 
magnetic Reynolds number is large in most parts of the com- 
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Fig. 6. Side by side logarithmic density contours and magnetic 
field lines (left panel) for the two opposite propagating jets, 
i.e. NRP (on the left) and NRA (on the right). Contours of V z 
and streamlines (right panel) for the same models. Model NRA 
displays higher vertical velocity around the transition region as 
compared to NRR 
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Fig. 7. Radial profile of log(p), V z and \og(pV z ) at z = 160 for 
models NIP (left; dotted line), NIA (left; dashed dotted), NRP 
(right; dotted) and NRA (right; dashed dotted). Resistive effects 
are negligible for the dynamics of the parallel case, while they 
amplify the features observed in the ideal anti-parallel one. 



putational domain, i.e. R m » 1, indicating that the diffusion of 
the magnetic field is not of particular importance. However, note 
that this definition of R m does not take into account the magnetic 
field reversal whose characteristic length is much smaller than r 
and hence the above argument excludes this region. Therefore, 
approximating resistivity with a small constant value is a useful 
approach to allow reconnection in our models without signifi- 
cantly affecting the validity of the ideal MHD results of the inner 
and outer regions. 

Moreover, we observe that the condition 10 < Rp < 100 
holds in most parts of our domain indicating that Ohmic heating 
is not important either. The value of Rp is also ill-defined in the 
anti-parallel case, since it is based on R m . Nevertheless, note that 
Rp is large at the axis where the plasma is hot and along the 
matching surface that the magnetic field reconnects, due to the 
high values of the plasma-/? there. 

Fig. [6] shows the distribution of log(p) and the magnetic field 
lines (left pair) of the final steady states of NRP (left) and NRA 
(right). The right pair displays the streamlines over the contours 
of V z for the same models. Initially, the asymmetric magnetic 
topologies lead to distinct current distributions. Reconnection is 
manifested at the regions where the magnetic field inverts sign, 
driving a different temporal evolution between the two hemi- 
spheres. Evidently, the presence of resistive effects in the anti- 
parallel case results in a spatial readjustment of the velocity pro- 
file in the inner regions. 

Another observed feature is the reversal of B z that takes place 
at the upper and outer domains of both NRP and NRA. This is 
a steady state configuration that originates from the resistive ef- 
fects at the base of the outflow and a shear in the velocity profile 
of the flow. In particular, although the two vector fields, B and 
V, are parallel at the bottom boundary, the low R m values sug- 



gests that diffusion dominates advection at the zones right above 
it. As a consequence, a small angle forms between the poloidal 
components of B and V before the magnetic Reynolds number 
of the flow acquires high values. Shearing effects increase this 
angle resulting in a negative B z component. 

Fig- El plots V z (middle) and the logarithms of p (top), and 
pV z (bottom) along the radial direction at z = 160 for mod- 
els NIP (left; dotted line), NIA (left; dashed dotted line), NRP 
(right; dotted line) and NRA (right; dashed dotted line). From 
the parallel setups we deduce that the prescribed value of resis- 
tivity has a minor effect on dynamics. On the contrary, non-ideal 
MHD effects amplify the features of the anti-parallel configura- 
tion, namely the decrease of the density at inner radii accompa- 
nied by an increase of the speed. In addition, the peak observed 
in p along the current sheet is more extended in NRA due to the 
physical reconnection applied. The mass flux profiles also reflect 
the asymmetry, despite the fact that p and V z change in opposite 
directions. Notably, beyond the matching surface all four models 
demonstrate similar behaviors. 

3.3. Time variable flows in parallel/anti-parallel configurations 

This section investigates the previously presented models when 
a time variable velocity is applied on the bottom boundary, lo- 
calized at the stellar component. The simulations are performed 
until the final time of the corresponding unperturbed cases is 
reached, during which several shocks propagate throughout the 
computational domain. 

Fig. [8] compares the logarithmic density contours and the 
magnetic field between models VRP (top) and VRA (bottom). 
The shock front is slightly faster in the anti-parallel configura- 
tion and seems to have a larger radius too. High matter concen- 
trations along the jet axis are observed at substantially different 
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Fig. 8. Logarithmic density and field lines for models VRP (top) and VRA (bottom). The two jets correspond to the two hemispheres 
of a YSO and are shown side by side for comparison. The shock front and matter condensations along the axis are found at different 
locations between the two cases. The stellar component of VRA displays a lower inner density as compared to VRP, whereas the 
inner part of its disk wind has higher values of p as matter accumulated along the current sheet is blown outwards from radial 
shocks. 
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Fig. 9. Logarithmic density profiles along the axis for the ideal 
MHD models (bottom), i.e. NIP (dotted line), NIA (dashed dot- 
ted), VIP (solid), VIA (dashed), and the resistive MHD models 
(top), NRP (dotted), NRA (dashed dotted), VRP (solid), VRA 
(dashed). The density has on average lower values in the anti- 
parallel cases and also peaks at different locations as compared 
to the parallel configuration. Especially in the presence of resis- 
tivity, decollimation drops the densities by an order of magnitude 
and significantly amplifies the differences between the parallel 
and anti-parallel models. 



heights, i.e. z = 160 for VRP and z = 210 for VRA. Furthermore, 
higher density is also found beyond the matching fieldline of 
VRA suggesting that the disk wind is affected as well. In fact, 
the matter trapped along the current sheet is pushed radially out 
by the inner shocks that travel across the jet flow. 



Fig. [9] plots the logarithm of p along the z axis for the ideal 
(bottom) and resistive (top) cases. Dotted/dashed-dotted lines 
are used for the non- variable parallel/anti-parallel configura- 
tions, and solid/dashed for the variable outflows, respectively. 
Density has lower overall values in the anti-parallel topology 
whereas flow fluctuations alter substantially its vertical profile. 
In particular, the high concentration regions are located at dif- 
ferent heights as compared to the parallel setup, but without 
demonstrating any systematic lagging. Both features seem to be- 
come stronger when physical resistivity is included. On the one 
hand this implies the non-negligible role of numerical diffusiv- 
ity (required for stability), and on the other, the fact that larger 
r] values could lead to significantly different results between the 
two hemispheres. Note that since emissivity is proportional to 
the density squared, it is expected that the above discrepancies 
could have a direct impact on observations. The pressure fol- 
lows approximately the behavior of p and hence it is not shown. 
As a result, the temperature has an almost flat profile for the 
ideal cases, whereas Ohmic heating slightly increases the aver- 
age temperature of the resistive models, although without any 
distinguishable effects among them. 

In Fig.[l0l which follows the notation of Fig.[9j the vertical 
velocity V z (top) is plotted versus the z direction. It was already 
anticipated from Fig. [7] that the parallel/anti-parallel configura- 
tions would show similar speeds. The small displacement of the 
shock front between the variable resistive models is evident here 
as well. However, having found discrepancies in the radial pro- 
file of V z , as well as in the density distribution, we introduce 
the following weighted average flow velocity, that incorporates 
emissivity in an approximate way: 



(V z ) 



j;-p 2 V z 2nrdr 



(10) 



This quantity provides a better estimate for the observational im- 
plications of the simulations than V z alone. The integration takes 
place within the cylinder of radius r max = 30, where p has its 
highest values. Note that for altitudes z > 200 the density peak 
along the current sheet is not inside the integrated domain. 

The bottom panels of Fig. [TO] display the vertical depen- 
dence of Eq. ([TO1) for the corresponding models. The ideal MHD 
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Fig. 10. Vertical velocity and weighted average of V z , Eq. (fTUb , 
along the z axis. The line style of each model follows from Fig.0 
Although the flow speed is the same among all cases, the quan- 
tity (V z ) does demonstrate asymmetries when non-ideal MHD 
effects are included. 

regime does not show any significant discrepancies beyond z - 
200. However, resistive cases do demonstrate velocity asymme- 
tries both in the variable and non- variable models. The difference 
is of the order of 30% and 20%, respectively, but since it depends 
on the specific 7] value assumed, it could be much higher or even 
negligible in a real astrophysical situation. In addition, note that 
including the current sheet within the integration domain results 
in a slower mean speed for the anti-parallel model. This can be 
seen by comparing the regions z < 200 and z > 200. In other 
words, the choice of r max also influences the absolute value of the 
estimated discrepancy and hence the resulting percentage should 
not be taken as a robust numbefl 

Although a parametric study of r] would clarify the asym- 
metry dependence on finite conductivity, such a task is beyond 
the scopes of this work. Besides, the actual resistivity value is 
an unknown factor in YSOs. It suffices that even in the case of 
7] = 1 .0, a small value that does not have significant effects on the 
parallel case, an asymmetry of -20-30% is found. Nevertheless, 
we have performed simulations with higher and lower resistivity 
values that seem to agree with the above conclusions. We also 
note that the vertical flow velocity is constant along z which in- 
dicates that there is no need to simulate larger spatial scales. 

3.4. External pressure 

All simulations of pressure-confined two-component jets reach 
quasi steady state configurations within a much larger simulation 
time as compared to the models presented above. This is due to 

3 On the contrary, our aim is qualitative and not quantitative, i.e. to 
show that velocity asymmetries can occur under some general consid- 
erations. A more complex average velocity could have been assumed, 
but we deliberately choose a simple function to constrain the number of 
unknown factors that could influence our results. 
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Fig. 11. Logarithmic density and field lines for models NIPL, 
VIPL, NIPH, VIPH, NIPVH and VIPVH (from left to right). The 
left and middle panels correspond to the northern and southern 
hemispheres of Fig. [2] respectively. The right panel displays the 
scenario of a very high external pressure in order to visualize 
the oscillation of the jet's cross section. Each panel shows the 
unperturbed model on the left and the time variable counterpart 
on the right. The outer structure as well as the condensations 
along the jet axis are static and unaffected by flow fluctuations. 



the fact that the initial conditions are farther away from the equi- 
librium described by the ADO and ASO solutions. Moreover, 
all final configurations possess static knots along the flow axis 
as well as a wave-like cylindrical outer surface that recollimates 
the flow at multiple locations. The separation of these high den- 
sity regions seems to increase the lower the ambient pressure 
applied. Notably, the variability of the inner flow does not seem 
to disrupt this structure, nor the shape or position of its surface. 
The jet models become asymptotically cylindrical, after hav- 
ing passed from a stage of oscillations in their radius, speed, 
density, and other physical parameters. Vlahakis & Tsinganos 
(119971) have shown that, under rather general assumptions, this 
is a common behavior of magnetized outflows which start non- 
cylindrically before they reach collimation. 

The logarithmic density along with the field lines are shown 
in Fig. HU for models NIPL, VIPL, NIPH, VIPH, NIPVH and 
VIPVH (from left to right). All of them have the same trun- 
cation radius, indicated in Fig. [3 but a different external pres- 
sure is imposed. Despite the pressure equilibrium holding right 
above the bottom boundary, the radial ram pressure of the flow 
expands the jet for small z- Eventually it becomes compensated 
and then recollimation occurs, compressing the flow at almost its 
initial diameter. At higher altitudes, the jet radius increases again 
and the process is repeated creating a sequence of static knots, 
a configuration that is stable and does not propagate with time. 
These features smooth out with z and the jet acquires a cylindri- 
cal shape far away, as seen in model NIPVH. We have verified 
the stability of the static knots with simulations carried out up 
to t = 500, though in a lower resolution. During time evolution, 
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Fig. 12. Vertical velocity contours for the models shown in 
Fig.CEU i.e. NIPL and VIPL (left panel), NIPH and VIPH (mid- 
dle panel), NIPVH and VIPVH (right panel). The velocity maps 
follow the density distributions. 



Kelvin-Helmholtz instabilities appear on the interface with the 
external medium, but they are gradually suppressed by the time 
the quasi steady state is reached. 

The jet structure of NIPH seems to be a smaller copy of 
NIPL, and the same holds true for NIPVH with respect to NIPH. 
This is due to the intrinsic self- similarity property of the ADO 
solution, namely, all field lines have the same shape but in a dif- 
ferent scale. Consequently, a similar outer structure is recovered 
in all pressure-confined cases, depending on the value of the 
external pressure imposed. Finally note that the mixing of the 
two analytical outflows introduces a length scale in the system 
and hence the above argument will not hold for arbitrarily small 
radii. 

Fig. [12] displays the velocity distribution for models NIPL 
(left pair; left), VIPL (left pair; right), NIPH (middle pair; left), 
VIPH (middle pair; right), NIPVH (right pair; left) and VIPVH 
(right pair; right). As expected, the velocity profile reflects the 
static oscillations of the jet's cross section since the mass flux is 
conserved. In addition, its radial dependence consists of layers of 
decreasing speed. Enforced time variability barely modifies the 
total distribution demonstrating the stability of the configuration. 
The small velocity observed in the outer medium is a product of 
the initial strong transient state in which the outflow is squeezed 
as the truncation surface (solid line, Fig. 0) moves upwards. 

The logarithmic density, the vertical velocity and (V z ) are 
plotted along the axis in Fig. [13] for the models NIPH (dotted), 
NIPL (dashed dotted), VIPH (solid) and VIPL (dashed). We no- 
tice that the discrepancies are stronger in the weighted veloc- 
ity, a quantity possibly relevant for observations. Moreover, (V z ) 
demonstrates a difference between the cases NIPH (dotted) and 
VIPH (solid) which is relatively smaller as compared to the vari- 
ations between the models NIPH and NIPL. This implies that 
the recollimation structure dominates, producing density peaks 
much higher than those of the enforced time variability. 
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Fig. 13. The quantities log(p), V z and (V z ) (Eq. [10| along the z 
axis, for models NIPH (dotted), NIPL (dashed dotted), VIPH 
(solid), VIPL (dashed). The pressure-confined jet simulations 
demonstrate velocity asymmetries in V z and become even more 
evident in the calculation of (V z ). 



Velocity asymmetries are evident in this class of simulations 
even without accounting emissivity. In particular, the middle plot 
of Fig. [T3l demonstrates a -10% deviation in V z between the two 
cases. The bottom panel suggests that the observed speed of the 
one jet could be twice the value of the other when weighted with 
the square of the density, at least locally. Note that in these cases 
the integration is calculated for r max = 20. Therefore, if a YSO 
outflow finds itself inside a non-uniform environment in which 
the pressure ratio between the two hemispheres is as low as 2, a 
significant velocity asymmetry could be measured. 

3.5. Mass loading 

A third possible scenario is that of an increased mass loading at 
the base of the outflow, such as a long lived asymmetric accretion 
or a strong radiation source heating unevenly the disk surface. 
In that case, the wind density can be different above and below 
the equator and so will the acceleration. However, on the one 
hand our numerical setup is not fully appropriate to investigate 
this mechanism in a consistent way, and on the other, the results 
of such simulations did not provide adequate evidence for the 
formation of velocity asymmetries. Therefore, we briefly discuss 
this case without presenting an extensive analysis. 

We have carried out simulations assuming the initial setup of 
NIP, but setting a two times higher density on the bottom bound- 
ary. These models allow us to examine the readjustment of the 
velocity profiles within a higher mass loading regime, naively 
assuming that all other jet properties remain the same at its base. 
Despite the significant modification of the boundary conditions, 
no considerable asymmetries have been found. Regarding the 
stellar component, we attribute this negative result to the fact that 
the outflow is already super- Alfvenic. A more consistent way 
to check how a high plasma mass will affect the velocity is to 
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take into account the acceleration regions. However, this is not a 
straightforward task since it is not clear how the energy input re- 
quired to drive the jet will change in a high mass loading regime. 
In addition, the much smaller characteristic time and lengths in- 
volved cannot be easily included in our comparably larger scale 
simulations. On the other hand, although the disk wind does in- 
clude sub-Alfvenic regions, they appear at large radii of low den- 
sity and speed that contribute negligibly to the velocity profile. 

To sum up, the large scale two-component jet models pre- 
sented in this paper, cannot fully capture the case of an asym- 
metric mass loading. The implicit assumption that the flow on 
the bottom boundary of each hemisphere has different mass but 
exactly the same velocity might not be valid if the accelera- 
tion mechanisms vary strongly with respect to the ejected mass. 
Nevertheless, we have verified that in the cases that this holds 
true, no substantial asymmetries are expected. 

4. Summary - conclusions 

In this paper we address the velocity asymmetries of YSO out- 
flows by investigating two classes of candidate mechanisms that 
could possibly generate this feature. 

The first class depends on the intrinsic properties of the YSO 
and assumes a parallel and an anti-parallel magnetic field con- 
figuration, one for each hemisphere. The application of physical 
reconnection results in different density and velocity distribu- 
tions between the two sides, leading to observable speed dis- 
crepancies. This scenario can provide asymmetric jets coming 
from intrinsic physical conditions and processes, even without 
considering the complex dynamics of the star-disk interaction. 
The relative questions that arise are whether multipolar fields in 
the star-disk system can exist and survive for a time scale com- 
parable to that of jet propagation, as well as what is the actual 
value of resistivity in YSO outflows, that could amplify or even 
suppress the asymmetry phenomenon. 

The second class of mechanisms is based on external effects, 
namely when the YSO resides in a non-uniform environment. 
Imposing distinct outer pressures at some boundary line of the 
jet is found to directly affect the degree of collimation of each 
flow, which in turn results in significantly modified propaga- 
tion speeds. This mechanism is appropriate to explain an out- 
flow asymmetry on large time scales. In addition, static knots are 
found to manifest in the jet structure due to multiple recolhma- 
tion locations. Their stability is robust despite the enforced flow 
perturbations, whereas their separation is found to be closely 
associated with the ambient pressure value. There is some ev- 
idence for stationary shocks in some systems (Matt & Bohm 
120031 Bonito et al. I20TT1 Schneider et al. |20TTT> , possibly asso- 
ciated with the collimation of the flow. However, shocks in pro- 
tostellar jets are generally observed to propagate with the flow 
(Reipurth & Bally 2001 and references therein). The fact that 
standing recollimation shocks, predicted in non- variable simu- 
lations with a confining pressure, are not generally observed in 
YSO jets implies that either the dynamical evolution of the out- 
flow is not strongly affected by the environment or that recolli- 
mation shocks are present but undetected. 

In the case of the RW Aur jet, it is not clear enough from the 
available observations whether the asymmetry is associated with 
the c entral engine (Woitas et al. 120021 Hartigan & Hillenbrand 
2009 ) or with environmental effects (Melnikov et al. l2009b . Both 
scenarios could in principle be feasible. We have arbitrarily se- 
lected the values of our parameters to explore both candidate 
mechanisms in a general manner without attempting to model 
this particular object. In order to understand the applicability of 



a particular model to interpret the velocity asymmetry, more de- 
tailed observations are required, as well as more advanced nu- 
merical models, including 3D geometry and radiative processes. 

Acknowledgements. We would like to thank an anonymous referee for help- 
ful suggestions and the A&A editor M. Walmsley for useful comments that 
lead to a better presentation of this work. We would also like to thank P. 
Tzeferacos, O. Te§ileanu and E. M. de Gouveia Dal Pino for fruitful discussions 
as well as J.-P. Chieze, F. Thais, C. Stehle and E. Audit for their support dur- 
ing the preparation of this paper. This research was supported by a Marie Curie 
European Reintegration Grant within the 7th European Community Framework 
Programme, ("TJ-CompTON", PERG05-GA-2009-249164). 



References 

Aburihan, M., Fiege, J., Henriksen, R., & Lery, T. 2001, MNRAS, 326, 1217 
Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 
Bessolaz, N., Zanni, C, Ferreira, J., Keppens, R., & Bouvier, J. 2008, A&A, 478, 
155 

Bonito, R., Orlando, S., Miceli, M., et al. 2011, ApJ, 737, 54 
Cemeljic, M., Gracia, J., Vlahakis, N., & Tsinganos, K. 2008, MNRAS, 389, 
1022 

Chagelishvili, G. D., Bodo, G., & Trussoni, E. 1996, A&A, 306, 329 
Coffey, D., Bacciotti, E, Woitas, J., Ray, T. P., & Eisloffel, J. 2004, ApJ, 604, 
758 

Daou, A. G., Johns-Krull, C. M., & Valenti, J. A. 2006, AJ, 131, 520 

Delamarter, G., Frank, A., & Hartmann, L. 2000, ApJ, 530, 923 

Donati, J.-E, & Landstreet, J.D. 2009, ARA&A, 47, 333 

Edwards, S., Fischer, W., Hillenbrand, L., & Kwan, J. 2006, ApJ, 646, 319 

Fendt, C. 2006, ApJ, 651, 272 

Fendt, C. 2009, ApJ, 692, 346 

Ferreira, J., & Casse, F. 2004, ApJ, 601, 139 

Ferreira, J., Dougados, C, & Cabrit, S. 2006, A&A, 453, 785 

Frank, A., & Noriega-Crespo, G. 1994, A&A, 290, 643 

Frank, A., & Mellema, G. 1996, ApJ, 472, 684 

Gracia, J., Vlahakis, N., & Tsinganos, K. 2006, MNRAS, 367, 201 

Hartigan, P., & Hillenbrand, L. 2009, ApJ, 705, 1388 

Konigl, A. 1982, ApJ, 261, 115 

Kwan, J., Edwards, S., & Fischer, W. 2007, ApJ, 657, 897 
Long, M., Romanova, M. M., & Lovelace, R. V. E. 2007, MNRAS, 374, 436 
Long, M., Romanova, M. M., & Lovelace, R. V. E. 2008, MNRAS, 386, 1274 
Lovelace, R. V. E., Romanova, M. M., Ustyugova, G. V., & Koldoba, A. V. 2010, 

MNRAS, 408, 2083 
Matsakos, T., Tsinganos, K., Vlahakis, N., et al. 2008, A&A, 477, 521 
Matsakos, T., Massaglia, S., Trussoni, E., et al. 2009, A&A, 502, 217 (M09) 
Matt, S., & Bohm K.-H. 2003, PASP, 115, 334 
Matt, S., Winglee, R., & Bohm, K.-H. 2003, MNRAS, 345, 660 
Matt, S., & Pudritz, R. 2005, ApJ, 632, 135 
Matt, S., & Pudritz, R. 2008, ApJ, 678, 1109 
Meliani, Z., Casse, E, & Sauty, C. 2006, A&A, 460, 1 

Melnikov, S. Yu, Eisloffel, J., Bacciotti, E, Woitas, J., & Ray, T. P. 2009, A&A, 
506, 763 

Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 
Mignone, A., Zanni, C, Tzeferacos, P., et al. 2012, ApJS, 198, 7 
Ouyed, R., Clarke, D., & Pudritz, R. E. 2003, ApJ, 582, 292 
Perin, M. D., & Graham, J. R. 2007, ApJ, 670, 499 
Reipurth, B., & Bally, J. 2001, ARA&A, 39, 403 
Salmeron, R., Konigl, A., & Wardle, M. 2011, MNRAS, 412, 1162 
Sauty, C, & Tsinganos, K. 1994, A&A, 287, 893 
Sauty, C, Trussoni, E., & Tsinganos, K. 2002, A&A, 389, 1068 
Sauty, C, Meliani, Z., Lima, J. J. G., et al. 2011, A&A, 533, 46 
Schneider, P. C, Gunther, H. M., & Schmitt, J. H. M. M. 2011, A&A, 530, 123 
Stute, M., Tsinganos, K., Vlahakis, N., Matsakos, T., & Gracia, J. 2008, A&A, 
491, 339 

Tesjleanu, O., Massaglia, S., Mignone, A., Bodo, G., & Bacciotti, F. 2009, A&A, 
507,581 

Tesjleanu, O., Mignone, A., Massaglia, S., & Bacciotti, F. 2012, ApJ, 746, 96 
Trussoni, E., Tsinganos, K., & Sauty, C. 1997, A&A, 325, 1099 
Tzeferacos, P., Ferrari, A., Mignone, A., et al. 2009, MNRAS, 400, 820 
Valenti, J. A., & Johns-Krull, C. M. 2004, Ap&SS, 292, 619 
Vlahakis, N., & Tsinganos, K. 1997, MNRAS, 292, 591 
Vlahakis, N., & Tsinganos, K. 1998, MNRAS, 298, 777 

Vlahakis, N., Tsinganos, K., Sauty, C, & Trussoni, E. 2000, MNRAS, 318, 417 
Wang, J. C. L., Sulkanen, M. E., & Lovelace R. V. E. 1992, ApJ, 390, 46 
Woitas, J., Ray, T., Bacciotti, E, Davis, C, & Eisloffel, J. 2002, ApJ, 580, 336 
Yang, H., Johns-Krull, C. M., & Valenti, J. A. 2007, AJ, 133, 73 



11 



T. Matsakos et al.: Velocity asymmetries in YSO jets 

Zanni, C, Ferrari, A., Rosner, R., Bodo, G., & Massaglia, S. 2007, A&A, 469, 
811 

Zanni, C, & Ferreira, J. 2009, A&A, 508, 1117 



12 



